Análisis de Datos de Área: El Analfabetismo en Nepal

Authors

Camilo Gómez

Sebastián Olarte

Descripción Del Conjunto De Datos: International Aid to Nepal (2007–2014)

El conjunto de datos utilizado corresponde al proyecto International Aid to Nepal (2007–14), desarrollado por el Center for Spatial Data Science (GeoDa Center). Su objetivo es integrar información geoespacial, socioeconómica y financiera que permita estudiar la pobreza, el desarrollo humano y la distribución de la ayuda internacional en los 75 distritos administrativos de Nepal. Este dataset combina datos espaciales y estadísticos.

Resumen General del Dataset

Característica Detalle
Título International Aid to Nepal (2007–14)
Observaciones 75 distritos (un polígono por distrito)
Variables 61 variables
Cobertura temporal 1991–2013 (ayuda internacional: 2007–2014)
Fuentes principales AidData y Open Nepal
Tipo de dato espacial Polígonos distritales (MULTIPOLYGON)
CRS WGS84 – EPSG:4326
Archivos shapefile .shp, .shx, .dbf, .prj, .qpj

Categorías y Descripción de Variables

Las 61 variables se organizan en cinco grupos temáticos principales.

A. Variables de Identificación Geográfica

Variable Descripción Unidad/Tipo
id Identificador único del distrito Numérica
name_1 Región administrativa Texto
name_2 Zona administrativa Texto
district Nombre del distrito Texto
lon Longitud del centroide Grados
lat Latitud del centroide Grados

B. Indicadores de Pobreza y Economía

Variable Descripción Unidad
DepEcProv Privación económica Índice
PovIndex Índice de Pobreza Humana Índice
PCInc Ingreso per cápita Rs.
PCIncPPP Ingreso per cápita ajustado por PPA USD PPP
PCIncMP Ingreso per cápita a precio de mercado Rs.

C. Indicadores Sociales, Salud y Servicios

Variable Descripción Unidad
Population Población total Personas
MalKids % de niños malnutridos (< 5 años) Porcentaje
Lif40 % que no alcanzará los 40 años Porcentaje
NoSafH20 % sin acceso a agua potable segura Porcentaje
VotNum Número de votantes registrados Personas

D. Indicadores Educativos (2012–2013)

Variable Descripción Unidad
BoyG1_5 Niños matriculados grado 1–5 Número
GirlG1_5 Niñas matriculadas grado 1–5 Número
KIDS1_5 Total niños grado 1–5 Número
SchoolCnt Número de escuelas Número
SCHLPKID Escuelas por 1.000 niños Ratio
SCHLPPOP Escuelas por 1.000 habitantes Ratio
AD_ILLIT Analfabetismo adulto (2011) Porcentaje
AD_ILGT50 1 si analfabetismo > 50%, 0 en otro caso Binaria

E. Variables de Ayuda Internacional

Cada sector contiene dos variables: xxCAMT (monto comprometido) y xxDAMT (monto desembolsado).

Código Sector Código Sector
AG Agricultura HEALT Salud
BANK Banca HUM Ayuda humanitaria
COMM Comunicaciones IND Industria
CON Conflicto MUL Multisectorial
EDU Educación SOC Infraestructura social
ENGY Energía TOUR Turismo
ENV Medio ambiente TRAN Transporte y almacenamiento
FOR Silvicultura WAT Agua y saneamiento
GOV Gobierno y sociedad civil TOT Total general

Lectura de Datos

Se realiza la lectura de los datos de Nepal.

Code
nepal <- st_read("nepal/nepal.shp", quiet = TRUE) # Se leen los datos de Nepal

nepal.st <- st_geometry(nepal) # Polígonos

#st_crs(nepal) # CRS

# Transformación a UTM

nepal.utm <- st_transform(nepal, crs = 32645) 

La Figure 1 muestra el mapa coroplético de la tasa de analfabetismo adulto en Nepal a nivel distrital. El mapa evidencia una marcada heterogeneidad espacial, con valores elevados concentrados principalmente en las regiones occidental y noroccidental del país, así como en parte de la franja sur, mientras que los niveles más bajos se observan en los distritos del centro y este.

La distribución espacial de los valores sugiere la presencia de agrupamientos de distritos con tasas similares de analfabetismo, indicando que la variable no se distribuye de manera aleatoria en el territorio. Este patrón visual constituye una evidencia descriptiva inicial que motiva la aplicación de métodos de análisis espacial para evaluar formalmente la dependencia espacial.

tmap_mode("view")

tm_shape(nepal.utm) +
  tm_polygons("ad_illit",
              fill.scale = tm_scale_intervals(n = 5, 
                                              style = "quantile", 
                                              values = "brewer.reds"),
              fill.legend = tm_legend(title = "Tasa de\nAnalfabetismo")) +
  tm_title("") +
  tm_layout(legend.position = c("left", "bottom"), frame = FALSE) +
  tm_scalebar(position = c("right", "bottom")) +
  tm_compass(position = c("right", "top"))
Figure 1: Mapa de la tasa de analfabetismo en Nepal por distrito

De esta forma, es posible ver que hay agrupaciones de distritos cuya tasa de analfabetismo es alta y otros para la cual la tasa es baja, indicando presencia de correlación espacial. Hacia el noroccidente y suroriente de Nepal se encuentran los distritos con una mayor tasa de analfabetismo.

Matriz de Vecindades

El primer paso consiste en definir los centroides de las áreas, que en este caso corresponden a los distritos de Nepal. Para ello, se emplean los centroides geométricos de cada distrito.

Code
coords <- st_coordinates(st_centroid(nepal.utm))

Pesos Espaciales

De esta forma, se procede con la definición de los vecinos para evaluar cuál es la matriz que mejor describe la relación espacial presente en los datos de analfabetismo de Nepal. En el caso de los vecinos físicos se obtienen las mismas conexiones para los métodos queen y rook.

Code
## Vecinos Físicos

nepal.nb <- poly2nb(nepal.utm, queen = TRUE)
nepal.lw_w <- nb2listw(nepal.nb, style = "W", zero.policy = TRUE)
nepal.lw_b <- nb2listw(nepal.nb, style = "B", zero.policy = TRUE)

## Vecinos basados en grafos

tri.nb <- tri2nb(coords)
tri.nb_w <- nb2listw(tri.nb, style = "W", zero.policy = TRUE)
tri.nb_b <- nb2listw(tri.nb, style = "B", zero.policy = TRUE)


soi.nb <- graph2nb(soi.graph(tri.nb, coords))
soi.nb_w <- nb2listw(soi.nb, style = "W", zero.policy = TRUE)
soi.nb_b <- nb2listw(soi.nb, style = "B", zero.policy = TRUE)


relative.nb <- graph2nb(relativeneigh(coords))
relative.nb_w <- nb2listw(relative.nb, style = "W", zero.policy = TRUE)
relative.nb_b <- nb2listw(relative.nb, style = "B", zero.policy = TRUE)


gabriel.nb <- graph2nb(gabrielneigh(coords), sym=TRUE)
gabriel.nb_w <- nb2listw(gabriel.nb, style = "W", zero.policy = TRUE)
gabriel.nb_b <- nb2listw(gabriel.nb, style = "B", zero.policy = TRUE)


## Vecinos basados en distancia

knn.nb_1 <- knn2nb(knearneigh(coords, k = 1))
knn.nb_2 <- knn2nb(knearneigh(coords, k = 2))
knn.nb_3 <- knn2nb(knearneigh(coords, k = 3))
knn.nb_4 <- knn2nb(knearneigh(coords, k = 4))


knn.nb_1_w <- nb2listw(knn.nb_1, style = "W", zero.policy = TRUE)
knn.nb_1_b <- nb2listw(knn.nb_1, style = "B", zero.policy = TRUE)


knn.nb_2_w <- nb2listw(knn.nb_2, style = "W", zero.policy = TRUE)
knn.nb_2_b <- nb2listw(knn.nb_2, style = "B", zero.policy = TRUE)


knn.nb_3_w <- nb2listw(knn.nb_3, style = "W", zero.policy = TRUE)
knn.nb_3_b <- nb2listw(knn.nb_3, style = "B", zero.policy = TRUE)


knn.nb_4_w <- nb2listw(knn.nb_4, style = "W", zero.policy = TRUE)
knn.nb_4_b <- nb2listw(knn.nb_4, style = "B", zero.policy = TRUE)

Luego de haber definido las matrices de conectividad, se realiza la elección de la mejor con ayuda del test I de Morán: se realiza el test y se escoge la que de un menor \(p\)-valor. Primero se almacenan todas las matrices en una lista y se define una función para realizar el procedimiento.

El sistema de hipótesis que se quiere juzgar es el siguiente, donde \(I\) representa el índice de Morán (Moraga (2023)):

\[ \begin{cases} H_0: I= \mathrm{E}[I] & \text{(no hay autocorrelación espacial)}\\ \text{vs} \\ H_1: I\not=\mathrm{E}[I] & \text{(hay autocorrelación espacial)}\\ \end{cases} \]

mat <- list(nepal.lw_w, nepal.lw_b,
            tri.nb_w, tri.nb_b,
            soi.nb_w, soi.nb_b,
            relative.nb_w, relative.nb_b,
            gabriel.nb_w, gabriel.nb_b,
            knn.nb_1_w, knn.nb_1_b,
            knn.nb_2_w, knn.nb_2_b,
            knn.nb_3_w, knn.nb_3_b,
            knn.nb_4_w, knn.nb_4_b)


names(mat) <- c(paste0("phys_", c("W","B")),
                paste0("tri_",  c("W","B")),
                paste0("soi_",  c("W","B")),
                paste0("rel_",  c("W","B")),
                paste0("gab_",  c("W","B")),
                paste0("knn1_", c("W","B")),
                paste0("knn2_", c("W","B")),
                paste0("knn3_", c("W","B")),
                paste0("knn4_", c("W","B")))

matrix_eval <- function(var, weights) {
  
  aux <- numeric(length(weights)) 
  
  for (i in seq_along(weights)) {
    aux[i] <- moran.mc(var,
                       weights[[i]], nsim = 9999,
                       alternative = "two.sided")$p.value
  }
  
  best_idx <- which.min(aux)
  
  list(
    best_index  = best_idx,
    best_weight = weights[[best_idx]],
    best_name = names(weights)[which.min(aux)],
    best_moran  = moran.mc(var, weights[[best_idx]], nsim = 9999,
                           alternative = "two.sided"),
    all_pvalues = aux
  )
}

Se ve cuál fue la matriz que mejor describe la correlación espacial de los datos.

Code
best.w <- matrix_eval(nepal.utm$ad_illit, mat)
best.weights <- best.w$best_weight

best.weights
Characteristics of weights list object:
Neighbour list object:
Number of regions: 75 
Number of nonzero links: 366 
Percentage nonzero weights: 6.506667 
Average number of links: 4.88 

Weights style: W 
Weights constants summary:
   n   nn S0       S1       S2
W 75 5625 75 32.85108 309.2524

Como se puede ver, esta corresponde a la obtenida considerando los cuatro vecinos más cercanos con el estilo W, es decir, con las filas estandarizadas. El mapa con esta matriz de conectividad se puede ver en la Figure 2.

nb_lines <- nb2lines(best.weights$neighbours, coords = coords, as_sf = TRUE)

st_crs(nb_lines) <- st_crs(nepal.utm)

ggplot() +
  geom_sf(data = nepal.utm, fill = "grey95", color = "grey80", size = 0.2) +
  # Conexiones
  geom_sf(data = nb_lines, color = "#2c3e50", size = 0.3, alpha = 0.6) +
  
  # Centroides
  geom_point(data = coords, aes(x = X, y = Y), 
             color = "#e74c3c", size = 1, alpha = 0.8) +
  
  # Estética
  theme_void() +
  labs(title = "Estructura de Conectividad Espacial",
       subtitle = paste0("Matriz de Pesos: ", best.w$best_name))
Figure 2: Mapa de Nepal con la Mejor Estructura de Conectividad Espacial

Se realizó el Test I de Moran mediante Monte Carlo con \(n=9999\) permutaciones. La distribución empírica del índice de Moran puede observarse en la Figure 3.

Figure 3: Histograma de los valores simulados del índice de Moran. La línea cortada representa del valor del índice de Moran observado.

De esta forma, se evidencia que al 5% de significancia se rechaza \(H_0: I = \mathrm{E}[I]\). Es decir, hay suficiente evidencia para concluir que hay presencia de correlación espacial con respecto a la tasa de adultos analfabetas en Nepal.

En la Figure 4 se muestra el diagrama de Moran de la tasa de adultos analfabetas en Nepal.

Figure 4: Dispersograma de Moran para la tasa de adultos analfabetas en Nepal
  • Índice de Moran (\(I = 0.595\)):
    El valor del estadístico es positivo y sustancialmente mayor que cero, lo que indica la presencia de una autocorrelación espacial global positiva fuerte*. Esto sugiere que los distritos geográficamente cercanos tienden a presentar tasas de analfabetismo similares: áreas con niveles altos tienden a agruparse espacialmente, al igual que áreas con niveles bajos.
    Dado que el \(p\)-valor es prácticamente cero (\(p < 0.001\)), se rechaza la hipótesis nula de ausencia de autocorrelación espacial. En consecuencia, el patrón observado no es producto del azar, sino que refleja una estructura espacial estadísticamente significativa.

    1. Cuadrantes

El diagrama se divide en cuatro cuadrantes definidos por las medias de la tasa de analfabetismo (eje horizontal) y su rezago espacial (eje vertical). Cada cuadrante representa un tipo particular de asociación espacial:

Cuadrante Patrón Descripción Ejemplos en la Figura Implicación Geográfica
Arriba–Derecha Alto–Alto (H–H) Distritos con tasas de analfabetismo altas rodeados por vecinos con valores igualmente altos. Casos en el extremo superior derecho. Hotspots de analfabetismo.
Abajo–Izquierda Bajo–Bajo (L–L) Distritos con tasas de analfabetismo bajas rodeados por vecinos con valores bajos. Lalitpur y puntos cercanos al origen. Coldspots de analfabetismo.
Arriba–Izquierda Bajo–Alto (L–H) Distritos con tasas bajas rodeados por vecinos con tasas altas. Surkhet. Outliers espaciales (rompen el patrón regional).
Abajo–Derecha Alto–Bajo (H–L) Distritos con tasas altas rodeados por vecinos con tasas bajas. Kapilbastu. Outliers espaciales con comportamiento diferenciado.

La pendiente positiva de la recta de regresión —equivalente al valor del Índice de Moran— y la concentración de observaciones en los cuadrantes Alto–Alto y Bajo–Bajo confirman que la tasa de analfabetismo adulto es un fenómeno espacialmente dependiente en Nepal. Este resultado sugiere que factores subyacentes al analfabetismo (por ejemplo, inversión en educación, pobreza estructural o condiciones socioculturales) no se distribuyen de manera aleatoria, sino que operan a escala regional y trascienden los límites administrativos distritales. En consecuencia, el diseño de políticas públicas orientadas a la reducción del analfabetismo debería adoptar un enfoque regional o interdistrital, concentrando esfuerzos en los hotspots identificados y analizando los outliers espaciales para comprender mecanismos locales y replicar experiencias exitosas.

Prueba de Autocorrelación Local

Con el fin de identificar clusters y evaluar la similitud entre áreas y sus vecinos de forma local, se realiza el test I de Moran de forma locaL. En la Figure 5 se muestran dos mapas en los que se observa autocorrelación espacial significativa en los distritos de las zonas noroccidental, central y suroriental de Nepal. En estas zonas predominan relaciones tipo alto–alto y bajo–bajo, lo que indica que los distritos con altas tasas de analfabetismo en adultos tienden a estar rodeados por distritos con tasas igualmente altas, y de igual forma para las tasas bajas.

lmoran <- localmoran(nepal.utm$ad_illit, 
                     best.weights, 
                     alternative = "two.sided")
nepal.utm$lmp <- lmoran[, "Pr(z != E(Ii))"]
Figure 5: Análisis de Autocorrelación Espacial Local para la tasa de analfabetismo.

Modelos de Regresión

Se ajustan distintos modelos de regresión, teniendo en cuenta las variables descritas al inicio del documento.

data <- as.data.frame(nepal.utm)

data <- data %>% 
  mutate(ag = AGDAMT / AGCAMT,
         engy = log((ENGYCAMT + 1 )/ (ENGYDAMT + 1)),
         gov = (GOVDAMT + 1) / (GOVCAMT + 1),
         healt = HEALTDAMT / population,
         mult = MULTDAMT / MULTCAMT,
         soc = (SOCDAMT + 0.001) / (SOCCAMT + 0.001),
         tot = TOTDAMT / population,
         log.schlppop = log(schlppop))

Modelo de Regresión MCO

Primero, se muestra el modelo de regresión que tuvo un mejor ajuste con respecto al \(R^2\) (0.4917) y la significancia de los coeficientes. Este modelo tiene la forma \[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \]

Code
form <- as.formula("ad_illit ~ pcinc + boyg1_5 + lif40 + soc + tot")

fit.OLS <- lm(form,
              data = data)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 57.9031 11.9913 4.8288 0.0000
pcinc -0.0231 0.0043 -5.4282 0.0000
boyg1_5 0.0002 0.0001 3.5416 0.0007
lif40 1.5015 0.5215 2.8790 0.0053
soc -25.6299 12.2164 -2.0980 0.0396
tot 0.0096 0.0049 1.9531 0.0549

Además, se realiza el Test I de Moran con los residuales del modelo, rechazando la ausencia de autocorrelación espacial al 5% de significancia. Es decir, los residuales presentan una dependencia espacial que no se pudo explicar con el modelo.

lmmorantest <- lm.morantest(fit.OLS, listw = best.weights)
t.lmmorantest <- tidy(lmmorantest)

kable(t.lmmorantest[,c("statistic", "p.value", "method", "alternative")])
statistic p.value method alternative
6.216516 0 Global Moran I for regression residuals greater

De esta forma, se hace necesario ajustar modelos de regresión espacial.

SAR - Spatial Autoregressive Model (Spatial Lag)

El modelo SAR incorpora la interdependencia espacial de la variable espacial con el \(\rho \in [-1,1]\). Este modelo tiene la estructura: \[ \mathbf{y} = \rho \mathbf{W} \mathbf{y} + \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \]

Code
fit.SAR <- lagsarlm(form,
                    data = data, listw = best.weights)

SEM - Spatial Error Model

Para el modelo SEM se tiene la estructura \[ \begin{aligned} \mathbf{y} =& \mathbf{X} \boldsymbol{\beta} + \mathbf{u}, \\ \mathbf{u} =& \lambda \mathbf{W} \mathbf{u} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \end{aligned} \] y este asume que hay factores no observados que generan la dependencia espacial. Es decir, la dependencia se incorpora en los errores.

Code
fit.SEM <- errorsarlm(form,
                      data = data, listw = best.weights)

SDM - Spatial Durbin Model

El modelo de Durbin integra tanto interdependencia espacial de la respuesta como efectos espaciales de las covariables. \[ \mathbf{y} = \rho \mathbf{W} \mathbf{y} + \mathbf{X} \boldsymbol{\beta} + \gamma \mathbf{W} \mathbf{X} + \boldsymbol{\varepsilon} , \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \]

Code
fit.Durbin <- lagsarlm(form,
                       data = data, listw = best.weights, type = "mixed")

SLX - Spatially Lagged X Model

En el modelo SLX se tienen en cuenta las relaciones entre las covariables de los vecinos. \[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \gamma \mathbf{W} \mathbf{X} + \boldsymbol{\varepsilon} , \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \]

Code
fit.D_lagX <- lmSLX(form,
                    data = data, listw = best.weights)

SDEM - Spatial Durbin Error Model

El modelo SDEM combina los modelos SEM y SLX. Este tiene la siguiente estructura \[ \begin{aligned} \mathbf{y} =& \mathbf{X} \boldsymbol{\beta} + \lambda \mathbf{W} \mathbf{X} + \mathbf{u}\\ \mathbf{u} =& \gamma\mathbf{W} \mathbf{u} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \end{aligned} \]

Code
fit.SDEM <- errorsarlm(form,
                       data = data, listw = best.weights, etype = "emixed")

Manski Model

Es un modelo complejo de la forma \[ \begin{aligned} \mathbf{y} =& \rho \mathbf{W} \mathbf{y} + \mathbf{X} \boldsymbol{\beta} + \mathbf{W} \mathbf{X} \theta + \mathbf{u}, \\ \mathbf{u} =& \lambda \mathbf{W} \mathbf{u} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \end{aligned} \]

Code
fit.Manski <- sacsarlm(form,
                       data = data, listw = best.weights, type="sacmixed")

SAC (SARAR)- Combined Spatial Autocorrelation Model

Por último, se tiene una combinación de modelos SAR y SEM. Es decir, se asume interdependencia espacial de la variable dependiente y factores no observados. \[ \begin{aligned} \mathbf{y} =& \rho \mathbf{W} \mathbf{y} + \mathbf{X} \boldsymbol{\beta} + \mathbf{u} \\ \mathbf{u} =& \lambda \mathbf{W} \mathbf{u}+ \boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon} \sim \mathrm{MVN}(\mathbf{0}, \sigma^2 \mathbf{I}) \end{aligned} \]

Code
fit.SARAR <- sacsarlm(form,
                      data = data, listw = best.weights, type="sac")

Selección del Modelo

Para realizar la selección del mejor modelo, se recurre a comparar los AIC.

Table 1: Comparación del AIC de los modelos ajustados
Modelo gl AIC
fit.SARAR 9 499.7997
fit.SDEM 13 501.6684
fit.Durbin 13 501.7228
fit.SEM 8 503.7218
fit.SAR 8 504.1058
fit.D_lagX 12 531.7039
fit.OLS 7 539.9456

En la Table 1 se puede observar que el modelo con un menor AIC corresponde al SARAR.

Interpretación

Debido a que el modelo SAC/SARAR fue el de menor AIC, se muestra su resúmen. Primero, se puede ver que no todas las variables son significativas, a diferencia del modelo de regresión ajustado por mínimos cuadrados ordinarios.

coef(summary(fit.SARAR)) %>% 
  kable(format = "html", digits = 4, align = 'c', escape = FALSE)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 32.4605 11.3319 2.8645 0.0042
pcinc -0.0183 0.0032 -5.6597 0.0000
boyg1_5 0.0001 0.0000 1.6876 0.0915
lif40 0.9806 0.4439 2.2090 0.0272
soc -10.1646 8.2507 -1.2320 0.2180
tot 0.0059 0.0036 1.6328 0.1025

Por otro lado, se presentan las estimaciones los parámetros \(\rho\) y \(\lambda\).

params.sac <- data.frame(Parámetro = c("$\\rho$", "$\\lambda$"),
                         Estimación = c(fit.SARAR$rho, fit.SARAR$lambda),
                         SE = c(fit.SARAR$rho.se, fit.SARAR$lambda.se))

params.sac %>% 
  kable(format = "pandoc", digits = 4, align = 'c', row.names = FALSE)
Parámetro Estimación SE
\(\rho\) 0.4682 0.1891
\(\lambda\) 0.5683 0.2027

De esta forma, con \(\hat{\rho}=0.468\) se confirma que hay dependencia espacial con respecto a la tasa de adultos analfabetas en Nepal. Más concretamente, un distrito puede afectar ``positivamente’’ al valor de sus vecinos, justo como se veía en la Figure 5. Por otro lado, \(\hat{\lambda} = 0.568\) indica que hay factores no incluídos en el modelo que presentan dependencia espacial. En este caso, ambas estimaciones son significativas al 5%.

Por lo tanto, se podría ver que a medida que los ingresos per cápita aumentan en un distrito, su tasa de analfabetismo disminuye. Por otro lado, si en un distrito el porcentaje de personas que no sobreviven hasta los 40 años aumenta, se espera que la tasa de analfabetismo aumente también.

Validación

Correlación Espacial

La primera evaluación que se realiza corresponde a ver si los residuales del mejor modelo ajustado siguen presentando autocorrelación espacial o no. Para esto, se realiza el Test I de Moran, concluyendo que, a una significancia del 5%, no se rechaza la hipótesis de ausencia de autocorrelación espacial de los residuales.

data$residuals <- resid(fit.SARAR)

moran.sac <- moran.test(data$residuals, listw = best.weights, alternative = "two.sided")
kable(tidy(moran.sac)[,c("statistic", "p.value", "method", "alternative")])
statistic p.value method alternative
0.2946557 0.7682569 Moran I test under randomisation two.sided

Homoscedasticidad

El siguiente supuesto que se prueba es el de varianza homogénea. Primero, en la Figure 6 se puede ver que los residuales no tienen ningún comportamiento que indique heterogeneidad de la varianza.

Figure 6: Residuales vs Valores Ajustados para el modelo SARAR

Además, se usa el test de Breusch-Pagan para juzgar \(H_0: \text{la varianza es homogénea}\).

bp.sac <- bptest.Sarlm(fit.SARAR)
kable(tidy(bp.sac)[,c("statistic", "p.value", "method")])
statistic p.value method
3.93796 0.5583824 studentized Breusch-Pagan test

Así, al 5% de significancia no se rechaza la hipótesis nula de homoscedasticidad.

Normalidad

En la Figure 7 se puede ver que los residuales tienden a estar cerca de la línea teórica, indicando un ajuste bueno.

Figure 7: QQ-plot de los residuales del modelo SARAR

Para corroborar esto, se realiza el test de Anderson–Darling, concluyendo que no se rechaza la hipótesis de normalidad de los residuales al 5% de significancia.

ad.sarar <- ad.test(data$residuals)
kable(tidy(ad.sarar))
statistic p.value method
0.6159227 0.105192 Anderson-Darling normality test

Estimación

A continuación, en la Figure 8 se presentan las estimaciones de la tasa de analfabetismo en adultos en Nepal.

nepal.utm$fitted <- fitted(fit.SARAR)
nepal.utm$residuals <- resid(fit.SARAR)
Figure 8: Mapa de predicciones de la tasa de analfabetismo en Nepal por distrito

Otros Modelos

Se decide comparar los resultados con otros modelos que no incorporan parámetros de asociación espacial como se vió en la sección anterior.

Moran Eigenvalues

Para realizar el uso de los vectores propios de la matriz \(\mathbf{M}\mathbf{W}\mathbf{M}\), con \(\mathbf{M} = \mathbf{I} - \mathbf{X}(\mathbf{X}^t\mathbf{X})^{-1}\mathbf{X}^t\), se usa la función SpatialFiltering() que toma el subconjunto de vectores propios que reducen la correlación espacial de los residuales en el modelo de regresión MCO.

Code
SF <- SpatialFiltering(form, data = data, nb = best.weights$neighbours, 
                       style = "W", verbose = FALSE, ExactEV = TRUE)

fit.SF <- lm(ad_illit ~ pcinc + boyg1_5 + lif40 + soc + tot + fitted(SF),
             data = data)

Por lo tanto, también se presenta el resumen de este modelo, que además tiene un \(R^2=0.8076\).

Estimate Std. Error t value Pr(>|t|)
(Intercept) 57.9031 7.1238 8.1281 0.0000
pcinc -0.0231 0.0025 -9.1371 0.0000
boyg1_5 0.0002 0.0000 5.9614 0.0000
lif40 1.5015 0.3098 4.8462 0.0000
soc -25.6299 7.2576 -3.5315 0.0008
tot 0.0096 0.0029 3.2876 0.0017
fitted(SF)vec3 -38.4032 4.9944 -7.6892 0.0000
fitted(SF)vec11 -21.0116 4.9944 -4.2070 0.0001
fitted(SF)vec15 21.8552 4.9944 4.3759 0.0001
fitted(SF)vec1 12.0327 4.9944 2.4092 0.0193
fitted(SF)vec4 10.6444 4.9944 2.1313 0.0375
fitted(SF)vec2 -9.9101 4.9944 -1.9842 0.0521
fitted(SF)vec10 -10.3387 4.9944 -2.0701 0.0431
fitted(SF)vec12 -9.5249 4.9944 -1.9071 0.0616
fitted(SF)vec21 -13.9456 4.9944 -2.7922 0.0071
fitted(SF)vec19 11.4086 4.9944 2.2843 0.0262
fitted(SF)vec8 8.1828 4.9944 1.6384 0.1069
fitted(SF)vec20 -9.6081 4.9944 -1.9238 0.0595
fitted(SF)vec9 -6.9303 4.9944 -1.3876 0.1708
SF.moran <- lm.morantest(fit.SF, best.weights, alternative = "two.sided")
t.SFmoran <- tidy(SF.moran)

kable(t.SFmoran[,c("statistic", "p.value", "method", "alternative")])
statistic p.value method alternative
-0.0167557 0.9866315 Global Moran I for regression residuals two.sided
Code
nepal.utm$residuals.SF <- resid(fit.SF)

SF.bptest <- bptest(fit.SF)

scatter.smooth(resid(fit.SF) ~ fitted(fit.SF),
               xlab = "Valores Ajustados", ylab = "Residuales")
legend("bottomright", 
       legend=c(paste0("Breusch-Pagan:", round(SF.bptest$statistic, 4)), 
                paste0(expression(p), "-valor: ", round(SF.bptest$p.value, 4))),
       cex=1,
       bg='salmon')

Code
SF.normal <- ad.test(resid(fit.SF))

qqnorm(resid(fit.SF),
       xlab = "Cuantiles Teóricos",
       ylab = "Cuantiles Muestrales", 
       main = "QQ-plot de los residuales")
qqline(resid(fit.SF))

legend("bottomright", 
       legend=c(paste0("Anderson-Darling:", round(SF.normal$statistic, 4)), 
                paste0(expression(p), "-valor: ", round(SF.normal$p.value, 4))),
       cex=1,
       bg='salmon')

De esta forma, es posible ver que, por ejemplo, los residuales ya no presentan correlación espacial y el supuesto de normalidad se cumple. Lo único que no se cumple es la homoscedasticidad.

GAM

Otra opción es usar los modelos GAM. Una ventaja de ellos, es que permite modelar covariables que no tengan una relación lineal con la variable respuesta. La forma de incluir la correlación espacial es mediante las coordenadas: s(X, Y), como se muestra en Bivand, Pebesma, and Gómez-Rubio (2008).

Code
library(mgcv)

data <- cbind(data, coords)

fit.gam <- gam(ad_illit ~ s(pcinc, bs = "cr") + boyg1_5
                          + s(lif40, bs = "cr") + s(soc, bs = "cr")
                          + s(tot, bs = "cr") + s(X, Y), data = data)

La única forma de interpretar los suavizamientos es con los efectos parciales de las covariables sobre la variable respueta. De la misma forma, se presenta su diagnóstivco.

(a) Efecto parcial de pcinc
(b) Efecto parcial de lif40
(c) Efecto parcial de soc
(d) Efecto parcial de tot
(e) Efecto espacial estimado
Figure 9: Efectos parciales sobre la tasa de analfabetismo
nepal.utm$residuals.gam <- resid(fit.gam)

gam.moran <-  moran.test(nepal.utm$residuals.gam, best.weights, alternative =  "two.sided")

tgam.moran <- tidy(gam.moran)

kable(tgam.moran[,c("statistic", "p.value", "method", "alternative")])
statistic p.value method alternative
-2.492309 0.0126916 Moran I test under randomisation two.sided
Code
gam.bptest <- bptest(fit.gam)

scatter.smooth(resid(fit.gam) ~ fitted(fit.gam),
               xlab = "Valores Ajustados", ylab = "Residuales")
legend("bottomright", 
       legend=c(paste0("Breusch-Pagan:", round(gam.bptest$statistic, 4)), 
                paste0(expression(p), "-valor: ", round(gam.bptest$p.value, 4))),
       cex=1,
       bg='salmon')

Code
gam.normal <- ad.test(resid(fit.gam))

qqnorm(resid(fit.gam),
       xlab = "Cuantiles Teóricos",
       ylab = "Cuantiles Muestrales", 
       main = "QQ-plot de los residuales")
qqline(resid(fit.gam))

legend("bottomright", 
       legend=c(paste0("Anderson-Darling:", round(gam.normal$statistic, 4)), 
                paste0(expression(p), "-valor: ", round(gam.normal$p.value, 4))),
       cex=1,
       bg='salmon')

El supuesto de no correlación espacial de los residuales se rechaza con una significancia del 5%. Con la introducción de las coordenadas no se alcanza a elimianar la dependencia espacial de los residuales.

Comparación Predicciones

Se presentan las predicciones para los tres modelos ajustados en la Figure 10.

nepal.utm$fitted.SF <- fitted(fit.SF)
nepal.utm$fitted.gam <- fitted(fit.gam)
Figure 10: Comparación de las de predicciones de la tasa de analfabetismo en Nepal por distrito

References

Bivand, Roger S., Edzer J. Pebesma, and Virgilio Gómez-Rubio. 2008. Applied Spatial Data Analysis with r. Use r! Springer. https://doi.org/10.1007/978-0-387-78171-6.
Bohórquez, Martha. 2024. “Estadística Espacial y Espacio-Temporal Para Campos Aleatorios Escalares y Funcionales: Notas de Clase.” Universidad Nacional de Colombia.
Center for Spatial Data Science (CSDS). 2017. “International Aid to Nepal (2007–14).” https://geodacenter.github.io/data-and-lab/nepal/.
Moraga, Paula. 2023. Spatial Statistics for Data Science: Theory and Practice with r. Chapman & Hall/CRC Data Science Series. Chapman & Hall/CRC.
Rüttenauer, Tobias. 2024. “Spatial Data Analysis.” https://arxiv.org/abs/2402.09895.